function [output] = nd_icu(ptstruct,whichgroup)
%ND_ICU	perform novelty detection in ICU patients
%	[output] = nd_icu(ptsturct,groups) 
%	
%
%	Inputs:     ptstruct    patient structure
%               groups      group of patients to investigate dynamics
%		
%		
%	Outputs:    output
%		
%		
%	
%	See also runND, NDtoy

%	References: N/A (First Paper)
%	
%	

%	Copyright 2013 MAF Pimentel
%	This program is free software: you can redistribute it and/or modify
%	it under the terms of the GNU General Public License as published by
%	the Free Software Foundation, either version 3 of the License, or
%	(at your option) any later version.
%	
%	This program is distributed in the hope that it will be useful,
%	but WITHOUT ANY WARRANTY; without even the implied warranty of
%	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%	GNU General Public License for more details.
%	
%	You should have received a copy of the GNU General Public License
%	along with this program.  If not, see <http://www.gnu.org/licenses/>.


%	$Author: MAF Pimentel$
%	$Revision: 1.0.0.0$
%	$Date: 05-Jun-2013 10:11:02$
%	Contact: marco.and.pimentel@gmail.com
%	Originally written on: PCWIN64

if nargin < 2
    whichgroup = 'mort28d';
end

if nargin < 1
    load '../data/ptstruct.mat';
end

LIMS = [20 350; ...             % HR
    0 300; ...                  % Pulse Pressure
    0 50; ...                   % Lactate
    0 1000];                    % Urine Output

PLOTLIMS = [40 250 1; ...       % HR
    0 200 1; ...                % Pulse Pressure
    0 50 1; ...                 % Lactate
    0 300 1];                   % Urine Output

LABELS = {'Heart Rate (bpm)','Pulse Pressure (mmHg)', ...
    'Acid Lactic','Urine Output (cc/kg/hr)'};

clc; close all;

DOSAVE = 0;
dirout = '../outputs/nd/';

%% Get ids for different classes of patients

[ids] = get_clids(ptstruct,whichgroup,1,0);

if isempty(ids)
    fprintf('\nCheck group of patients\n');
    return;
end

%% Cleaning data
%  for urine output and lactate
allUO = [];
for pt = 1 : length(ptstruct)
    if ~isempty(ptstruct(pt).data{1}.uo)
        allUO = [allUO; ptstruct(pt).data{1}.uo(:,2)];
    end  
end

allLA = [];
for pt = 1 : length(ptstruct)
    if ~isempty(ptstruct(pt).data{1}.lactate)
        allLA = [allLA; ptstruct(pt).data{1}.lactate(:,2)];
    end  
end


%% Conduct univariate comparisons distributions from different points 
%% of the patient trajectory
% Get useful data first

time = 12; % las hours and time to get data from
DATA = get_datalast(ptstruct,ids,time*60);
COLORG = {'g','r'};
for v = 1 : 4
    figure;
    for c = 1 : 2
        data = DATA{c}(:,v);
        data(isnan(data)) = [];
        data(data < LIMS(v,1)) = [];
        data(data > LIMS(v,2)) = [];
        subplot(2,1,c);
        hist(data,PLOTLIMS(v,1):PLOTLIMS(v,3):PLOTLIMS(v,2));
        h = findobj(gca,'Type','patch');
        set(h,'FaceColor',COLORG{c},'EdgeColor','k');
        xlabel(LABELS{v}); ylabel('Counts');
        fprintf('\nGroup %d, %s = %.1f (%.1f)',c,LABELS{v},mean(data),std(data));
        xlim([PLOTLIMS(v,1) PLOTLIMS(v,2)]);
    end
    disp('');
end

fprintf('\n');

%% Construct two-dimensional model of normality
% Select paramaters and filter training data
params = [1 2];

% Randomly select training and testing patients in the normal group
Nnor = length(ids{1});
Nabn = length(ids{2});

[iTrain, iTest] = crossvalind('HoldOut',Nnor, 0.3);
trainpts = ids{1}(iTrain); trainpts = trainpts(:)';
testnorpts = ids{1}(iTest); testnorpts = testnorpts(:)';
testabnpts = ids{2}; testabnpts = testabnpts(:)';

data = DATA{1}(iTrain,params);
[r,~] = find(isnan(data)); data(unique(r),:) = [];
indrm1 = find(data(:,1) < LIMS(1,1)); indrm2 = find(data(:,1) > LIMS(1,2));
indrm3 = find(data(:,2) < LIMS(2,1)); indrm4 = find(data(:,2) > LIMS(2,2));
indrm = [indrm1(:); indrm2; indrm3; indrm4];
data(indrm,:) = [];

% Perform normalization of training data
meanval = mean(data); stdval = std(data);
datan = data - repmat(meanval,size(data,1),1);
datan = datan./repmat(stdval,size(data,1),1);

% Construct model
addpath('netlab/');
MIX = Train_kde(datan);

% Analyse raw timeseries or 'normalisation'
dotimenorm = 0;

if ~dotimenorm
    % Get test data
    t = 4*60;                   % Define length of windows
    lastday = 1:13;
    nwin = t:t:lastday(end-1)*24*60;        % Define number of windows
    idslastday{lastday(end)} = [];
    
    % Split timeseries
    datatrain{length(trainpts)} = [];
    datatestnor{length(testnorpts)} = [];
    datatestabn{length(testabnpts)} = [];
    
    for i = 1 : length(trainpts)
        for w = nwin
            dataprov = get_datawin(ptstruct,trainpts(i),[w-nwin(1) w]);
            datatrain{i} = [datatrain{i}; w dataprov(params) -1];
        end
        % LoS
        if floor(ptstruct(trainpts(i)).los) == 0
            idslastday{1} = [idslastday{1}; i 1];
        elseif floor(ptstruct(trainpts(i)).los) > 12
            idslastday{13} = [idslastday{13}; i 1];
        else
            d = floor(ptstruct(trainpts(i)).los);
            idslastday{d} = [idslastday{d}; i 1];
        end
            
        % remove missing values
        [r,~] = find(isnan(datatrain{i}));
        datatrain{i}(unique(r),:) = [];
        % determine likelihood
        datatrain{i}(:,2:end-1) = datatrain{i}(:,2:end-1) - repmat(meanval,size(datatrain{i},1),1);
        datatrain{i}(:,2:end-1) = datatrain{i}(:,2:end-1)./repmat(stdval,size(datatrain{i},1),1);
        datatrain{i}(:,end) = Out_kde(datatrain{i}(:,2:end-1),MIX);
    end
    
    for i = 1 : length(testnorpts)
        for w = nwin
            dataprov = get_datawin(ptstruct,testnorpts(i),[w-nwin(1) w]);
            datatestnor{i} = [datatestnor{i}; w dataprov(params) -1];
        end
        % LoS
        if floor(ptstruct(testnorpts(i)).los) == 0
            idslastday{1} = [idslastday{1}; i 2];
        elseif floor(ptstruct(testnorpts(i)).los) > 12
            idslastday{13} = [idslastday{13}; i 2];
        else
            d = floor(ptstruct(testnorpts(i)).los);
            idslastday{d} = [idslastday{d}; i 2];
        end
        
        % remove missing values
        [r,~] = find(isnan(datatestnor{i}));
        datatestnor{i}(unique(r),:) = [];
        % determine likelihood
        datatestnor{i}(:,2:end-1) = datatestnor{i}(:,2:end-1) - repmat(meanval,size(datatestnor{i},1),1);
        datatestnor{i}(:,2:end-1) = datatestnor{i}(:,2:end-1)./repmat(stdval,size(datatestnor{i},1),1);
        datatestnor{i}(:,end) = Out_kde(datatestnor{i}(:,2:end-1),MIX);
    end
    
    for i = 1 : length(testabnpts)
        for w = nwin
            dataprov = get_datawin(ptstruct,testabnpts(i),[w-nwin(1) w]);
            datatestabn{i} = [datatestabn{i}; w dataprov(params) -1];
        end
        % LoS
        if floor(ptstruct(testabnpts(i)).los) == 0
            idslastday{1} = [idslastday{1}; i 3];
        elseif floor(ptstruct(testabnpts(i)).los) > 12
            idslastday{13} = [idslastday{13}; i 3];
        else
            d = floor(ptstruct(testabnpts(i)).los);
            idslastday{d} = [idslastday{d}; i 3];
        end
        
        % remove missing values
        [r,~] = find(isnan(datatestabn{i}));
        datatestabn{i}(unique(r),:) = [];
        % determine likelihood
        datatestabn{i}(:,2:end-1) = datatestabn{i}(:,2:end-1) - repmat(meanval,size(datatestabn{i},1),1);
        datatestabn{i}(:,2:end-1) = datatestabn{i}(:,2:end-1)./repmat(stdval,size(datatestabn{i},1),1);
        datatestabn{i}(:,end) = Out_kde(datatestabn{i}(:,2:end-1),MIX);
    end
    
    % Determine average values for model of normality
    % Compute average in which time window?
    
    output.averagedresults = zeros(length(nwin),7);
    output.averagedresults_1 = zeros(length(nwin),27);
    output.averagedresults_2 = zeros(length(nwin),27);
    output.averagedresults_3 = zeros(length(nwin),27);
    for w = 1 : length(nwin)
        output.averagedresults(w,1) = nwin(w)/(60*24);
        
        dataprov = [];
        for pt = 1 : length(datatrain)
            ind = find(datatrain{pt}(:,1) == nwin(w));
            if ~isempty(ind)
                dataprov = [dataprov; -log(datatrain{pt}(ind,end))];
            end
        end
        output.averagedresults(w,2) = mean(dataprov);
        output.averagedresults(w,3) = std(dataprov)/sqrt(length(dataprov));
        
        dataprov = [];
        for pt = 1 : length(datatestnor)
            ind = find(datatestnor{pt}(:,1) == nwin(w));
            if ~isempty(ind)
                dataprov = [dataprov; -log(datatestnor{pt}(ind,end))];
            end
        end
        output.averagedresults(w,4) = mean(dataprov);
        output.averagedresults(w,5) = std(dataprov)/sqrt(length(dataprov));
        
        dataprov = [];
        for pt = 1 : length(datatestabn)
            ind = find(datatestabn{pt}(:,1) == nwin(w));
            if ~isempty(ind)
                dataprov = [dataprov; -log(datatestabn{pt}(ind,end))];
            end
        end
        output.averagedresults(w,6) = mean(dataprov);
        output.averagedresults(w,7) = std(dataprov)/sqrt(length(dataprov));
        
        
        output.averagedresults_1(w,27) = nwin(w)/(60*24);
        for d = 1 : 13
            ptList = idslastday{d}(idslastday{d}(:,2) == 1,1);
            dataprov = [];
            for pt = ptList(:)'
                ind = find(datatrain{pt}(:,1) == nwin(w));
                if ~isempty(ind)
                    dataprov = [dataprov; -log(datatrain{pt}(ind,end))];
                end
            end
            output.averagedresults_1(w,2*d-1) = mean(dataprov);
            output.averagedresults_1(w,2*d) = std(dataprov)/sqrt(length(dataprov));
        end
        
        output.averagedresults_2(w,27) = nwin(w)/(60*24);
        for d = 1 : 13
            ptList = idslastday{d}(idslastday{d}(:,2) == 2,1);
            dataprov = [];
            for pt = ptList(:)'
                ind = find(datatestnor{pt}(:,1) == nwin(w));
                if ~isempty(ind)
                    dataprov = [dataprov; -log(datatestnor{pt}(ind,end))];
                end
            end
            output.averagedresults_2(w,2*d-1) = mean(dataprov);
            output.averagedresults_2(w,2*d) = std(dataprov)/sqrt(length(dataprov));
        end
        
        output.averagedresults_3(w,27) = nwin(w)/(60*24);
        for d = 1 : 13
            ptList = idslastday{d}(idslastday{d}(:,2) == 3,1);
            dataprov = [];
            for pt = ptList(:)'
                ind = find(datatestabn{pt}(:,1) == nwin(w));
                if ~isempty(ind)
                    dataprov = [dataprov; -log(datatestabn{pt}(ind,end))];
                end
            end
            output.averagedresults_3(w,2*d-1) = mean(dataprov);
            output.averagedresults_3(w,2*d) = std(dataprov)/sqrt(length(dataprov));
        end
        
    end
    
    figure; hold on
    errorbar(output.averagedresults(:,1),output.averagedresults(:,2), ...
        output.averagedresults(:,3),'Color',[0 1 0]);
    errorbar(output.averagedresults(:,1),output.averagedresults(:,4), ...
        output.averagedresults(:,5),'Color',[0 .5 0]);
    errorbar(output.averagedresults(:,1),output.averagedresults(:,6), ...
        output.averagedresults(:,7),'Color',[.5 0 0]);
    hold off
    legend('Train','Normal Test','Abnormal Test'); xlabel('Time (mins)'); ylabel('Novelty Score');
    
%     for d = 1 : 13
%         figure; hold on;
%         errorbar(output.averagedresults_1(:,27),output.averagedresults_1(:,2*d-1), ...
%             output.averagedresults_1(:,2*d),'Color',[0 1 0]);
%         errorbar(output.averagedresults_2(:,27),output.averagedresults_2(:,2*d-1), ...
%             output.averagedresults_2(:,2*d),'Color',[0 .5 0]);
%         errorbar(output.averagedresults_3(:,27),output.averagedresults_3(:,2*d-1), ...
%             output.averagedresults_3(:,2*d),'Color',[.5 0 0]);
%         hold off
%         legend('Train','Normal Test','Abnormal Test'); xlabel('Time (mins)'); ylabel('Novelty Score');
%         pause;
%     end
else    
%     disp('');
end

% Get test data for display
% Make a movie
if length(params) == 2
    DOANIM1 = 1;
    figure;
    [xgrid, xgrida, xgridb] = get_grid2D(datan, 100, 100);
    ygrid_vec = Out_kde(xgrid, MIX);
    ygrid = reshape(ygrid_vec, 100, 100);
    surface(xgrida, xgridb, ygrid,'EdgeColor','interp'); axis tight; zlim([0 max(ygrid(:))+.1]); box on;
    view([70 40]);
    if DOANIM1
        hold on
        [C,h] = contour3(xgrida, xgridb, ygrid, [0.01 max(ygrid(:))+.1]);
        set(h,'LineWidth',5,'CData',[0 0 0],'EdgeColor','k');
        thresh = .01:.01/10:max(ygrid(:))-.05;
        thresh = [thresh fliplr(thresh)];
        for i = 1 : length(thresh)
            set(h,'Visible','off');
            [~,h] = contour3(xgrida, xgridb, ygrid, [thresh(i) max(ygrid(:))+.1]);
            set(h,'LineWidth',3,'CData',[0 0 0],'EdgeColor','r');
            drawnow;
            pause(.02);
        end
    end
    hold off
end

if length(params) == 2
    DOANIM2 = 1;
    figure;
    [xgrid, xgrida, xgridb] = get_grid2D(datan, 100, 100);
    ygrid_vec = Out_kde(xgrid, MIX);
    ygrid = reshape(ygrid_vec, 100, 100);
    surface(xgrida, xgridb, ygrid,'EdgeColor','interp'); axis tight; zlim([0 max(ygrid(:))+.1]); box on;
    view([-90 90]);
    if DOANIM2
        hold on
        % select an abnormal patient
        timeseries1 = datatestnor{2}(:,2:4);
        
        % select a normal patient
        timeseries2 = datatestabn{4}(:,2:4);
        
        for i = 2 : size(timeseries1,1)
            plot3(timeseries1(1:i,1),timeseries1(1:i,2), ...
                timeseries1(1:i,3),'-og','MarkerSize',2,'LineWidth',2);
            drawnow;
            pause(.5);
        end
        
        for i = 2 : size(timeseries2,1)
            plot3(timeseries2(1:i,1),timeseries2(1:i,2), ...
                timeseries2(1:i,3),'-ok','MarkerSize',2,'LineWidth',2);
            drawnow;
            pause(.5);
        end
    end
    hold off
end



disp('');

end